function [SOL, d] = get_equivalentBEDtreatment(param, delta_t, free, t_f1, t_f2, D_BED_0, IT, n, LQL, vd, use_Markov)
    
% INFO: Function to evaluate different RT treatments with the same BED

%% OUT:

    %SOL: (vec) evolution of tumor volumes for each treatment
    %d: dose per fraction for each treatment with the equivalent BED
    
%% IN:

    % param: (vec) parameters
    % delta_t: (real) time step
    % free: (vec) controls free growth and start of treatment condition
    %    free(1): (int) 1 for free behavior until an specific time or volume, 0 otherwise 
    %    free(2): (int) 1 for time, 0 for volume
    %    free(3): (real) initial volume or time
    % t_f1: (real) maximum time allowed to reach the initial state
    % t_f2: (real) simulation time after initial state is reached
    % D_BED_0: (vec) [n fractions, d dose per fraction]
    % IT: =0 no anti-CTLA4, otherwise IT
    % n: (vec) number of fractions for each treatment
    % LQL: =1 if LQL, otherwise (modified) LQ (with PARAM(33)) 
    % vd: =1 if vascular death (D > 15) 
    % use_Markov: =1 if Markov, otherwise deterministic
    
%% Algorithm

%Calculate BED_0
BED_0 = 1/(param(3)/param(4)) * D_BED_0(1) * D_BED_0(2)^2 + D_BED_0(1) * D_BED_0(2);
m = length(n);

%Calculate the dose per fraction for each treatment given by the number of fractions and the equivalent BED
d = zeros(1, m);
for i = 1:m
    d(i) = (-n(i) + sqrt(n(i)^2 - 4 * n(i)/(param(3)/param(4)) * (-BED_0))) / (2 * n(i) / (param(3)/param(4)));
end

SOL = zeros(m, round(t_f2/delta_t)+1);
tum_cells = zeros(m, round(t_f2/delta_t)+1);

for j = 1:m

    %Dose vector
    D = d(j) * ones(1, n(j));
    t_rad = 0 : (length(D) - 1);
    t_treat_c4 = [2 5 8];
    if IT == 0
        param(25) = 0;
        t_treat_c4 = 0;
    end
    param(35) = 0;
    t_treat_p1 = 0;

    [sol, t_eq, ~, ~, N] = radioimmuno_response_model(param, delta_t, free, t_f1, t_f2, D, t_rad, t_treat_c4, t_treat_p1, LQL, vd, use_Markov);
    SOL(j,:) = sol(round(t_eq/delta_t) : round(t_eq/delta_t) + round(t_f2/delta_t));
    tum_cells(j,:) = N(round(t_eq/delta_t) : round(t_eq/delta_t) + round(t_f2/delta_t));
   
end